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CN ; Abstract 

A dynamical system description of the transition process in sliear flows with no linear instability 
^O ■ starts with a knowledge of exact coherent solutions, among them travelling waves (TWs) and 

relative periodic orbits (RPOs). We describe a numerical method to find such solutions in pipe flow 
<^' and apply it in the vicinity of a Hopf bifurcation from a TW which looks to be especially relevant 

h _ 

^H , modulated streaks. This RPO, like the TW from which it bifurcates, sits on the laminar-turbulent 

O ' boundary separating initial conditions which lead to turbulence from those which immediately 

J>^' relaminarise. 
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for transition. The dominant structural feature of the RPO solution is the presence of weakly 
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I. INTRODUCTION 

Ever since the experiments of Osborne Reynolds in 1883,- it has been known that the 
flow of a Newtonian fluid in a circular pipe undergoes an abrupt transition from a laminar to 
a turbulent state. This transition is governed by one dimensionless parameter, the Reynolds 
number Re := Ud/u {U is the bulk velocity, d the pipe diameter and u the kinematic 
viscosity of the fluid). Transition is usually observed experimentally for Re ~ 2000 ± 250 
but can be delayed to larger values in especially well-controlled experiments.-i^ The laminar 
state (Hagen-Poiseuille flow) is believed linearly stable for all Re, so that transition to 
turbulence cannot be explained by an ordered sequence of bifurcations starting from the basic 
state. Instead, pipe flow belongs to a more general class of parallel shear flows, including 
plane Couette and plane Poiseuille flow, which undergo flnite-amplitude instability.- The 
transitional dynamics of these flows then appears organised around the presence of other 
solutions disconnected from the laminar state.-'^'^'^ Such new solutions are expected to have 
a relatively simple structure compared to the turbulent dynamics and to be unstable. These 
expectations are conceptually satisfying for two reasons. Firstly, relatively simple three- 
dimensional mechanisms (e.g. the 'self-sustaining process') have been identified in low-Re 
flows^*^ allowing turbulence to be sustained over long times, which evokes the possibility 
of organised dynamics. Secondly, instability is needed to rationalise why neither steady nor 
periodic exact states are ever found as limiting states to emerge from an initially turbulent 
flow. The fact that statistically recurrent states have been clearly observed, e.g. streaky 
velocity flelds with an intrinsic wavelength in most wall-bounded shear flows, ^ is most easily 
explained by imagining repeated visits in phase space to the vicinity of exact solutions of 
the governing equations. 

Progress in numerical computation has flrst led to the discovery of some exact solutions 
in plane shear flows^ii^ii^ and more recently pipe flow.— ^^^ They are all equilibria (steady 
states) or relative equilibria (travelling waves which are steady in a moving frame) and all 
possess a small number of unstable eigendirections. They appear through saddle-node bifur- 
cations and are disconnected from the laminar state. There is a tremendous interest in these 
solutions as very similar structures have been observed transiently in experiments.— i^^ It is 
worth emphasizing that their simple dynamics (steadiness or constant speed propagation) is 
inherent to the method used to seek them, and undoubtedly hides a larger variety of other 



solutions. At some point the physical description of the flow in terms of three-dimensional 
structures ceases to be enlightening and it is useful to think in terms of a trajectory in phase 
space. Formally, the Navier-Stokes equations can be projected onto any inflnite-dimensional 
basis of incompressible velocity flelds to yield an autonomous dynamical system of the form 
dX/dt = F{X,Re). The choice of the basis, and hence the particular projection, is not 
crucial. The topology of the phase space associated to pipe flow near transition is thought 
to be organised around one locally attractive point (the 'laminar' flow) and a chaotic sad- 
dle, i.e. a set of unstable solutions and their heteroclinic and homoclinic connections.-i^^ 
From this point of view, it is necessary to know which states the saddle is made of in order 
to understand the dynamics near transition. It has been demonstrated that the unstable 
solutions known so far in pipe flow, the travelling waves (TWs), are visited only for about 
0(10-20%) of the time.— li^ This indicates that, provided the picture is correct, more solu- 
tions of a different type have to be sought. The next level up in the hierarchy of solutions 
from relative equilibria is relative periodic orbits (RPOs). They are an extension of periodic 
orbits (in the same way as traveling waves correspond to degenerate flxed points) due to the 
invariance of the equations with respect to the azimuthal and axial directions. Physically, a 
RPO corresponds to a flow which repeats itself periodically in a given moving frame, trans- 
lating and/or rotating along the pipe axis at a constant rate. In the cylindrical coordinate 
system {s,9,z) aligned with the pipe axis, it is a time- dependent velocity fleld v{s,9,z,t) 
satisfying 

v{s, e,z,t + T) = v{s, 9 + Ae,z + Az, t) (1) 

for some constants T, A^ and Az. The initial motivation behind the search for periodic 
behaviour is Periodic Orbit Theory,— which states that any dynamical average of a smooth 
chaotic dynamical system can be evaluated in a deterministic way by flnding all unstable 
periodic orbits of the system up to a certain period, and carefully averaging over them. 
Consideration of continuous symmetries in the original PDEs - here translation in z and 
rotation in 6' - extends this expansion over all relative periodic orbits. ^^ Even if complete 
knowledge of all periodic orbits of the system seems computationally ambitious today, a 
start needs to be made in developing numerical tools to flnd them in anticipation of greater 
computing power tomorrow. Some progress has already been made in plane Couette flow 
with the recent discovery of periodio^^ and relative periodic orbits^i embedded in turbulence. 



albeit in small flow domains. Beyond the obvious computational savings, restricting the 
largest wavelength in the flow has proved a very useful tool for isolating key dynamics^. 
In the same spirit, we will here concentrate on relatively short 'periodic' pipes only a few 
diameters long. 

In this study we are interested in RPOs which are located in the laminar-turbulent 
boundary, the invariant subset of phase-space on which trajectories neither relaminarise 
nor become turbulent. States belonging to this subset can be visited transiently by 
phase-space trajectories during the transition process. Phase-space trajectories constrained 
to remain precisely on this laminar-turbulent boundary have been found to approach a 
chaotic attractor™ centered on one particular TW solution, the 'asymmetric TW found by 
Pringle & Kerswell.— The importance of this particular TW has been confirmed by further 
calculations which show that trajectories constrained to lie on this laminar-turbulent 
boundary recurrently approach this TW.^^ We therefore concentrate on looking for RPOs 
which bifurcate from this asymmetric TW. 

The paper is organised as follows. Section [III explains in detail the numerical method 
chosen to find RPOs from a starting guess and section III describes the new branch of RPOs 
found. Subsection IIVI confirms that the RPO is on the laminar-turbulent boundary and then 
explores the likelihood of connections between the RPO and other TWs. Section IVl discusses 
the relevance of the RPO to the transition process as well as the numerical limitations of 
the method. 

II. NUMERICAL PROCEDURE 

A. Governing equations 

We consider the incompressible flow of Newtonian fluid in a cylindrical pipe and adopt 
the usual set of cylindrical coordinates (s, 6', z) and velocity components u = us + v0 -\- wz. 
Units of length and flow speed are taken as the pipe radius and the bulk speed U. The 
computational domain is (s, 6, z) G [0, 1] x [0, 27r] x [0, L], where L = 27i/a is the length of the 
pipe. The non-dimensionalised incompressible three-dimensional Navier-Stokes equations 



read 

— + {u-V)u = -Vp + -^^^u, (2) 

V ■ n = 0, (3) 

where Re := Ud/v. The flow is driven by a constant mass-flux condition 

p\ /'27r 

I/tt / / wsdsde = l, (4) 

Jo Jo 

as in recent experiments. -i2ii2^ Although time is calculated in units of d/2U, all times are 
quoted hereafter in the usual units of d/U. The boundary conditions are periodicity across 
the pipe length u{s, 9, z) = u{s, 9, z + L) and no-slip on the walls u{l, 9, z) = 0. In the 
non-dimensionalisation used, the expression of the Hagen-Poiseuille flow is uhp = 2(1 — s^)i. 



B. Time-stepping code 

The basic tool for the numerical determination of periodic orbits is an accurate time- 
stepping code. The direct numerical simulation code written by A. P. Willis^^'^iiS^ ^vas 
adopted which is based upon the toroidal-poloidal representation of the velocity field^ 

n = V X (\E'z) + V X V X ($£). (5) 

The two scalar potentials $ and \E' are discretised using high-order finite differences in the 
radial direction s, and with spectral Fourier expansions in the azimuthal 9 and axial z 
directions. At any discrete radial location Sj, {j = 1, ..., A^), each of the potentials (e.g. here 
$) is therefore discretised according to the formula: 

K M 

^s,,9,z,t-a)= Y, Yl ^,Ut)e'^'^'^''''^ (6) 

k=-K m=-M 

The resolution of a given calculation is described by a triplet {N,M,K). The phase-space 
associated to the Navier-Stokes equations is the set of complex coefficients {^jkm, ^jfcm}, 
which corresponds to a dynamical system with n ~ 8MNK real degrees of freedom (typically 
n = O(IO^)). Its metric is defined by the Euclidean scalar product (, ). A shift back in 
physical space by {Az, A9) corresponds in phase-space to the transformation: 



Time discretisation is of second-order, using a Crank-Nicholson scheme for the diffusion 
term and an Euler predictor step for the non-hnear terms. In this study we confine our 
calculations to a pipe of length L = 27r/0.75 ~ 8.38 radii = 4.19 c?. This length has been 
chosen because it is a wavelength for which the asymmetric TW^^ is known to be on the 
laminar-turbulent boundary and for which most data were already available. 



C. The Newton-Krylov method 



1. Periodic orbits 



A periodic orbit (of period T) of a dynamical system X = F{X) is sought as a solution 
(X,T) of the equation G{X) = 0, where 

G{X) := 0^(X) - X (8) 

Here (t>^{X) refers to the point at time t > on the trajectory starting from X at time 
t = 0. A value of the period t = T has to be chosen to uniquely define G{X). We choose 
T as the time minimising the JC— dependent scalar function g : t ^ |0*(X) — Xp over a 
given time interval. This is done by accurate interpolation of the zeros of its derivative, the 
scalar function g' : t ^ 2(0*(X) — X,d(f)'^{X)/dt). Figure [1] shows schematically that the 
time t = T is picked up along the trajectory when its tangent vector is orthogonal to the 
difference vector (f)^{X) — X, which ensures that the trajectory has come back closest to 
the initial point. A standard method to search for zeros of G in low-dimensional systems 
(typically n < O(IO^)) is the Newton-Raphson method. This iterative algorithm, based on 
successive Taylor expansions of G, produces a sequence of iterates X^''\ k >0 defined by 

X('=+i) = XW+(5XW, (9) 

where 

J{X^^^)dX^^^ = -G(X(^)), (10) 

J being the Jacobian matrix associated to G. When no analytical expression for J is 
available, it is generally computed by finite differences, at the cost of one evaluation of G 




FIG. 1: Schematic view of phase-space. Definition of the functional G to minimise in order to find 
exact periodic orbits. 

for each column of the matrix. In very high dimension n > O(IO^), however, the time 
needed for n + 1 evaluations of G, as well as the storage of J (not to mention solving 
equation fITOl) ). can become prohibitive, thus we turn our attention towards matrix-free 
Inexact Newton-Krylov methods. 

At every Newton step k, we solve the system (ITOl) using the GMRES algorithm.— This 
iterative linear solver only requires matrix- vector products that can be evaluated using finite 
differences, since for a vector p and e sufficiently small (e.g. 10^^) we have: 

G(X('=) + ep)-G(XW; 



Jp 



+ 0(e) 



(11) 



Hence no matrix needs to be stored explicitly. Using (ITT!) . GMRES iteratively builds an 
orthogonal basis of the Krylov subspace spanned by {rg, Jtq, J^Tq, ...}, starting from Vq = 
—G{X^^'). At each GMRES step, a Gram-Schmidt procedure reorthonormalises the Krylov 
subspace producing a basis {Vq, Vl, V2, ...}. Then an approximation to the solution of (TTOjl 
is constructed on this basis, until convergence is decided according to the criterion: 



|jW5X('=) + G(^)|<r7('=)|G('=)|. 



(12) 



The 'forcing' constant •r]'^'^'^ appearing in formula ( IT2l) is ideally zero, and numerically it 
should be in principle very small. However, non-vanishing values of 77 can allow for a much 



quicker convergence of the overall Newton scheme, by avoiding useless oversolving. Values of 
rj ~ O(10~^) have produced faster convergence and even better performances were observed 
when updating i]^^^ with the choice 2 of Eisenstat & Walker.— 

2. Double dogleg step 

Newton-Raphson algorithms are known to converge only if the initial guess is sufficiently 
close to a zero of the function G, which in high dimension often means no convergence at 
all. Moreover, even close to a solution, Newton-Raphson steps either can be too large or 
nearly orthogonal to the gradient of |Gp, in which case the algorithm stagnates and the 
classical strategy of linear backtracking (also called 'damped Newton') is of no help.™ To 
remedy this, it is useful to embed the algorithm into what is commonly called a 'globally' 
convergent strategy. Here 'global' does not mean that convergence is achieved whatever the 
starting point, but instead that the basin of attraction of a solution is reasonably enlarged. 
We adopt the 'double dogleg step' method proposed by Dennis and Schnabel,— which is a 
subclass of the 'trust region' algorithms (see Viswanath^ for an implementation of another 
subclass call the 'hook step' approach). All these methods start from the knowledge of the 
Newton step SX^ (here we omit the superscript (k) ), and aim at constructing a new vector 
6X whose norm is bounded by a 'trust length' 6c- The way dX is chosen is described in the 
next paragraph. At the first stage of the loop, 6X = 5Xn and 6c = {6X^1. Then at every 
iteration of the trust region algorithm, a check is made to see whether \G{X^'^' + dX)\ is 
sufficiently smaller than \G{X^''^)\. If it is, then the loop is finished and we set X^''^^^ to 
X^''^ + dX; if not, then 6c is divided by 2, a new vector SX is computed and the process is 
repeated. 

In the case of the double dogleg step, two descent directions for |Gp are used to generate 
the new vector dX : the first one is the Newton step SX^, the second is the steepest descent 
direction (see figure [2]): 

- ^V|G|2 = -J*G, (13) 

where J* is the transpose of the Jacobian matrix J. We can define the Cauchy Point 
X^''^ + SXqp, where the vector dXcp minimizes the quadratic form SX -^ \JSX + Gp 
along the steepest descent direction (whereas the overall minimizer of this quadratic form 
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defines the Newton step 5Xn)- Projection of SXcp onto tlie Krylov space (already built to 
find 6Xn at the same Newton iteration) allows a low-cost approximation of 5Xcp, as long 
as the dimension of the Krylov subspace itself is low. If the Vi s form an orthonormal basis 
of the Krylov subspace and K is its dimension, we have: 

K K 

V|G|2 ~ J2^J'G, Vi)V, = J2{G, JV^Vi. (14) 

i=l 4=1 

In (11^ . the scalar coefficients {G,JVi) are already known up for 1 < i < K since 
(normalized by — |^o|) they form the first row of the Hessenberg matrix used for the 
GMRES inversion.^ 

Finally, given dX^ and SXqp, the double dogleg algorithm finds the intersection between 
the piecewise linear curve linking the Cauchy Point to X + fidX^, and the ball of radius 
6c. So the new trial vector is chosen as: 

6X = SXcP + Xc (i^SXn - SXcp) . (15) 

where /i is a constant set to 0.8 (following ref.'^^) and Ac > is chosen such that |<5X| = 6c- If 
\6c\ < \5Xcp\, the dogleg curve along which optimisation is achieved is simply the steepest- 
descent direction of IGp. 

3. Relative periodic orbits 

Allowing for shifts in 6 and z implies a modification of the function G. If we seek a 
relative periodic orbit satisfying v{r, 6,z,t + T) = v{r,6 + A6, z + Az, t) for any t, we have 
to change G{X) into G{X) = 0!!^Ae-Az(^) ~ -^- ^^ this expression, 0-^(X) is shifted 
back in the 6- and z-directions using ([7]), before being compared to the initial point of the 
orbit. T is chosen as the time minimising the new |Gp, whereas A6 and Az are considered 
as two extra variables of the Newton algorithm:— we define a (n + 2)-dimensional vector 
X+ = (X, A9, Az) and the function G+ : 7?"+^ ^ ^n+2 g^^j^ ^j^^^^ 

G+(X+) = G,(X), (2 = l,..,n), (16) 

while the two last components of G~^ are defined by the scalar products in K^: 

G.|„(X-) = (G,|g\, (17) 




FIG. 2: Schematic view of the double dogleg technique (adapted from Dennis & Schnabel^^) in 
phase space. X^ > is the state at the k Newton iteration, 6Xj\f represents the Newton step 
computed using the inexact GMRES method, and — V|Gp is the steepest descent direction of 
the residual from X'^'^'. At a given step, a new state X^^^^> = JC^"^ + 6X is suggested as the 
intersection between the double dogleg curve (medium thick line) and the ball of radius 5c- 



Gn+2i^ 



G, 



dAz 



(18) 



These two derivatives are easily computed by central finite differences. Starting from a 
guess Xq (which includes a starting point in phase space as well as two initial values of A^ 
and Az), the Newton-Krylov method is now applied to locate zeros of the function G^. It 
is worth reiterating that although T looks to be on the same footing as A^ and Az, this 
is not the case. The condition < G, dcf)^ /dt >= which defines T is imposed at every 
Newton/dogleg iteration whereas the equivalent conditions for A6 and Az are only truly 
reached when convergence to a RPO has been achieved. 

Non-rotating TW solutions are a special case of RPOs, when T and Az are linked by the 
relation Az = cT, with c the axial propagation speed of the wave. Because of this degeneracy, 
TW solutions with zero azimuthal speed can be sought, using the same algorithm, by fixing 
the value of Az to 2n/a, and setting G^^2 = 0- 
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III. RESULTS 

A. Hopf bifurcation of traveling w^aves 

The Newton-Krylov algorithm described in Section [TTl was used to find a RPO bifurcating 
off the asymmetric TW branch in a periodic pipe of length L ~ 4.19 rf (a = 0.75) where it 
is known to be embedded in the laminar-turbulent boundary. The asymmetric TW branch 
originates from a supercritical bifurcation at Re ~ 1770 off a "mirror-symmetric" branch 
of TWs.~ The asymmetric TW propagates axially with a phase speed c = 1.34 f/ but does 
not rotate. For Re ~ 1785, it displays a pair of high-speed streaks sandwiching an eccentric 
low-speed streak (see Figure [3]). A streamwise vortex tube is present in the vicinity of the 
low-speed streak, ensuring the regeneration of the streak via a lift-up effect. The two high- 
speed regions extend to the other side of the cross-section, and another low-speed streak 
of lesser intensity is present on the right of the original one. As Re increases, the solution 
becomes more asymmetric and the streaks on the right of the Figure tend to disappear. 

The numerical resolution chosen in this work is (A^, M, K) = (86, 16, 5) to match that used 
by Pringle & Kerswell in the 6- and z-directions. This is sufficient to observe a drop-off of 6 
decades in the axial spectrum and more than 8 decades in the radial and azimuthal energy 
spectra. The rapid drop-off in the axial wavenumbers is typical of lower-branch solutions 
and their weakly wiggling streak structure. It is even more pronounced as Re — > cxo.— 
The Navier-Stokes equations have been linearised around the TW solution expressed in 
the Galilean frame moving with speed c along the axis. The corresponding eigenvalue 
problem was solved numerically using an Arnoldi routine for Re up to 5000.— There are 
4 unstable eigenvalues (2 real and 2 complex conjugate ones) for Re < Ren = 1785.6 and 
only 2 real ones for Re > Ren indicating a Hopf bifurcation at Re = 1785.6 (see Figure 
H]). The center manifold theorem states that for Re close enough to Ren, a periodic orbit 



must exist in the same moving frame. Its amplitude scales locally like 0{^^/\ReH — Re\).— 
The period of the orbit is given by the imaginary part of the Hopf eigenvalue pair, here 
T = 2tx / Iin{\H) ~ 43 d/U . A good estimation for the shifts is the distance travelled by the 
TW in the axial and azimuthal directions during this period, hence Az^°) = cT r^ 58.7 d and 
A9^^^ = (since this wave has no azimuthal phase speed). The starting point used for the 
Newton-Krylov algorithm is a slight perturbation of the TW along one of the two conjugate 
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FIG. 3: Velocity profile of the asymmetric TW solution at the Hopf bifurcation (a = 0.75, Re = 
1785). From left to right, from top to bottom : z/L = 0,0.25,0.5,0.75. The arrows indicate the 
cross-stream velocity, and the shading indicates the difference between the streamwise velocity and 
the laminar profile (light /white indicating positive values and dark/red negative values). 



neutral directions, denoted as the eigenvector en'. 



(19) 



For values of e ~ O(10~^) and Ren — 1 < Re < Ren, the algorithm converged to a RPO, 
within a residual r^o = \X_az(T) — X(0)|/|X(0)| of 0(10^^). The converged values of the 
period T and the shift Az are very close to the expected values, while the RPO appears 
not to rotate {A6 = 0). The new solution branch bifurcates towards the regime where the 
original branch has fewer stable directions, namely Re < Ren- This Hopf bifurcation is 
hence supercritical, although the new solution branch extends to lower values of Re. From 
now on we will also use the reduced positive parameter 6 = Ren — Re. 
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FIG. 4: Stability analysis of the asymmetric TW for a = 0.75 near the Hopf bifurcation point at 
Re = Ren = 1785.6: growth rates (in units of U/D) against Re. Each line either indicates the 
locus of a real eigenvalue (solid) or a complex conjugate pair (dashed) as Re changes. The blue 
(dark) lines correspond to those which are shift-&:-reflect symmetriol^ while the green (light) lines 
are anti-shift-&-reflect symmetric. The TW itself bifurcates off the mirror-symmetric TW branch 
near Re = 1770. 

B. Continuation of the RPO 

Continuation of these orbits along the i?e-axis is used to produce a bifurcation diagram. 
Once one RPO is known at a given value of 6, it is straightforward to use it as an initial 
condition for the Newton-Krylov algorithm at slightly different values of 6. Progression 
towards larger 6 (i.e. lower Re) by using small enough steps enables one in principle to 
track the solution down until the branch possibly bends back. 

Trajectories corresponding to RPO solutions have been traced on a two-dimensional plane 
{D', E') for various values of 5. Here D' := Re~^ / 1^ >< u'\'^<fx and E' := \ j \u'\'^d?x are 
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FIG. 5: Numerically found RPOs in the {D',E') plane for different values of (5 = 4 (solid red) 
5 = 22 (dotted blue) 6 = 47 (dashed green), with 5 = Ren — Re- The asymmetric (TWA) and 
mirror-symmetric (TWM) TW branches are also indicated (black dotted curves), parametrised by 
Re. The TWA at Re = Ren is indicated by the red arrow. The RPO corresponding to 5 = 47 
does not look perfectly closed, because of numerical issues detailed in Subsection IIIICI Note that 
the TWA branch exists only for Re > 1770 {5 < 15). 



respectively the dissipation and the kinetic energy of the disturbance u' := u — uhp to the 
laminar flow (see Figure [5]). Such a projection makes these RPOs look periodic, whereas in 
the non-moving frame they are only relative periodic orbits. The orbits are slightly elliptic 
and correspond to a slow modulation of the TW from which they have bifurcated, whereas 
a TW at a given Re would appear as a dot. 
We define their normalized amplitude by : 



5E 



Er, 



Er, 



K 



max ~r -t-'min 



(20) 



where the subscripts min and max denote extrema along a cycle of the RPO {5E = for 
a TW solution). The further 6 is away from the bifurcation point, the larger the amplitude 
of the RPO, as is expected from a regular Hopf bifurcation. The values of T and Az both 
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FIG. 6: Normalised amplitude {SE)"^ for the RPO as a function of Re (lower blue data), along with 
asymptote (SE)'^ = 2 x W^S where 6 = Ren — Re (blue dashed line). The minimal normalised 
residual r^o = \X^/^ziT) — X{0)\/\X{0)\ achieved by the Newton-Krylov algorithm is also shown 
as a function of Re (upper red data). Continuation was stopped when the numerical residual 
reached 10~^. 



keep the same order of magnitude. For S up to 50, the tendency for the orbits is to achieve 
excursions towards regions of higher dissipation and energy. The energy variations along one 
cycle stay small for the range of Re considered. Figure M shows the dependence of the orbit 
amplitude 6E on Re: 6E grows monotonically with the reduced parameter 6 and is well 
approximated by the expression 6E = 0(62) near the bifurcation point. The curve has been 
continued numerically until 6 = 47. Above this value, convergence becomes questionable - 
see Figure [6] - as discussed in Subsection IIIICI and in the conclusion. 

The value of Az indicates the distance travelled by the exact coherent structure as it 
propagates down the pipe in time T. This defines a phase speed for the RPO, 



crpo := Az/T, 



(21) 
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FIG. 7: Axial phase velocity c (in units of f7) vs. Re for the RPO (red), the asymmetric TW 
(TWA) and the mirror-symmetric TW (TWM) (both in black). 

corresponding to the axial speed of the frame in which the orbit would be exactly periodic. 
Figure [7] is a diagram showing the phase velocity of both the asymmetric TW branch and 
the mirror-symmetric TW branch, compared to the phase velocity of the RPO which varies 
little over the range of Re considered. Note that the RPO exists at values of Re below that 
at which the asymmetric TW appears (see Figure [T]). 

Figure [8] shows the evolution of the velocity field in a cross-section moving with the speed 
crpo at (5 = 47 (i.e. Re ~ 1738). This moving frame is chosen to emphasize slow-time 
variations. Near the left of the cross-section, the flow is strongly reminiscent of the pattern 
of the original asymmetric TW, which would give 4 identical slices. There is a slight periodic 
modulation of the shape and position of the low-speed streak. The high-speed streaks close 
to the wall look more robust in shape during a cycle, but their intensity fluctuates. As b 
increases, the intensity inside the patches of axial velocity on the right of Figure [8] fluctuates 
more in time than the ones of the left. The larger 5, the more symmetric the cross-sectional 
pattern looks (this feature is also shared with the asymmetric TW branch as it approaches 
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FIG. 8: Velocity field of the RPO in the moving frame, at a cross-section defined by z — cptpot = 
(see Eq. EH), for (5 = 47 {Re = 1738). From left to right, from top to bottom : t = 7.2, t = 18.7, 
t = 30.3 and t = 41.8 in units of d/U. The colours and arrows are as in Figure [3l 

the mirror-symmetric branch). The cross-section velocity is more steady than the axial 
velocity, except in the vicinity of the weaker low-speed streak on the right of the pipe, where 
the vortex looks 'attached' to the streak. 

Figure [9] displays the friction factor A, defined using dimensional variables as 



^-'lA^""^ 



(22) 



of the RPO branch. For any given value of Re, A along a cycle is always below those of the 
asymmetric and mirror-symmetric TW solutions. The highest value Amax reached during a 
cycle gets closer to that of the mirror-symmetric branch as Re is decreased. 
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FIG. 9: Friction factor A vs. Re for the RPO (red), the asymmetric TW and the mirror-symmetric 
TW (both black). In the case of the RPO, only the extremal values Amin and Amax over a cycle 
are shown. 

C. Instability of the RPO and numerical issues 

While performing the continuation in S, we observed that Tod (the value around which 
the numerical residual stagnates) increases steadily with 6. This trend is clear from Figure 
El While animation of the velocity field does not show any discontinuity in time for 6 > 25, 
the projection on the {D',E') plane no longer looks perfectly closed (see Figure E]). We 
deliberately stopped the continuation near 5 ~ 50 when it was difficult to reduce Too below 
10-2. 

In order to understand this behaviour of Voo, we first need an idea of how accurate the 
search for RPOs is expected to be. The minimum residual Too is limited by both the accuracy 
of the time-stepping scheme and the natural instability of the trajectory. Any component 
of any point on the numerically found RPO is accurate to at least 0{em), the machine 
precision. After a period T, this numerical discrepancy has increased by a factor e^^ , with 
A being now the largest Floquet exponent of the RPO. Indeed a useful rule-of-thumb for 
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estimating the lowest residual attainable is: 

O (ne„e^^) . (23) 



Tr. 



We can illustrate this formula by simply time-stepping from a sufficiently converged asym- 
metric TW at 5 = 0. After the TW has travelled one single pipe length, r(t = L/c) ~ 10~^°; 
whereas after a time T (the period of the RPO), r{T) ~ 10~^. Using expression (l23l) . we 
can estimate \tw ~ ln(10^)/T ~ 0.26f//(i. Sufficiently close to the Hopf bifurcation, both 
the RPO and the TW are expected to have the same leading Floquet exponent A, hence the 
estimation \rpo{^ = 0) ~ 0.26 U /d. Since the period T does not vary much (less than 1%) 
within the range of 5 considered, we can extract from the data in Figure 6 a trend for the 
Floquet exponent of the RPO. We found the following scaling : \rpo{,5) ~ 0{5'^), with 7 
very close to |. Unsurprisingly, the RPO is more unstable away from the bifurcation point. 
A possible explanation is the stronger slow-time modulation of the streaks (whose strong 
transversal velocity gradients are already responsible for the inertial instability of the TW) 
as 5 increases. This harmonic modulation of the streaky field over a full cycle is likely to 
increase the instability of the whole flow. 



IV. CONNECTION WITH A TW SOLUTION 

A leading question is how the RPO found in this paper is related to transition to turbu- 
lence in a pipe. The asymmetric TW is already known to play a special role. Firstly, it lies 
on the laminar-turbulent boundary (also called the 'edge of chaos— 1^), which is the set of 
initial conditions separating trajectories which directly relaminarise from those which lead to 
turbulent behaviour. This means that inflnitesimal perturbations to the TW can either lead 
to transition or relaminarisation, depending on their exact form. Secondly, recent studies 
have shown that for analogous parameters, but larger Re, this asymmetric TW is generally 
the only exact travelling wave solution closely visited by trajectories when constrained on 
the laminar-turbulent boundary. ^^ It is thus interesting to see if other kind of exact recurrent 
states, like the RPO found here fulfll the expectation of: (a) lying on the laminar-turbulent 
boundary as well; (b) connecting to other states via heteroclinic connections; and (c) getting 
transiently approached by transitional trajectories. 
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We have chosen to concentrate on S = 2 {Re = 1783.6), and perturbed the RPO by two 
random vectors in the following manner : 

X(t = 0) = Xrpo + e {coscj) ei + sincj) 62) . (24) 

In this expression, e is a small parameter (here 10"'^), ei and 62 are two random vectors 
of norm unity, which are generic enough to possess non-zero components along the unstable 
manifold of the RPO. is a shooting angle which has been varied in the range [0, 2tt]. We 
notice that several trajectories starting from the new state X{t = 0) relaminarise directly 
while other trajectories undergo a short 'turbulent' transient, indicated by a dramatic 
rise in the kinetic energy and dissipation of the disturbance to the laminar flow, and by 
a loss of symmetry of the flow. Hence the RPO found for Re = 1783.6 is also on the 
laminar-turbulent boundary. 

Establishing heteroclinic connections between the RPO and other exact coherent struc- 
tures is an involved and technically demanding task. We can, however, seek suggestive evi- 
dence for such links by exploring the evolution of trajectories originating close to the RPO 
and confined to remain in the laminar-turbulent boundary. This was attempted for the 
RPO at 5 = 2 by refining the angle defined in f l24|) to find a phase-space trajectory which 
stays near the laminar-turbulent boundary, i.e. neither directly relaminarises nor undergoes 
a turbulent transient. This is done via shooting method based on bisection of the value of 0. 
Refining up to three significant digits results in a trajectory called Hi, which eventually 
relaminarises after a significantly long transient along the laminar-turbulent boundary. For 
a short duration during the transient (less than 10 d/U), all velocity components oscillate 
with the same apparent period on a short-time scale, while the energy reaches a plateau. 
Based upon previous experienced^, we recognised this as the signature of a travelling wave 
solution lying nearby in phase space. The scalar function 

r^,„(t)= mini r,(t'), 1^ = 0} (25) 

where 

was calculated along the whole trajectory. This function, closely linked to |G|, measures 
how recurrent a flow is at a given location in space (ignoring the possibility of shifted 
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recurrences). Figure [TO] shows rmm as a function of the starting point on the edge trajectory 
Hi. While t < 60 d/U, Tmin is very low (0(10^^)) and the flow is nearly recurrent. 
The slow-time modulation reflects the modulation of the flow along the RPO. Later the 
trajectory leaves the neighbourhood of the RPO because of its instability, and rmin increases 
up to large values of ~ 0.2. At a later phase corresponding to the plateau in energy, Tminif) 
displays a clear dip (labelled Nl in Figure [TO!) down to values of ~ 5 x 10^^, and then 
increases again. Such low values of Tmin are never reached if the angle is chosen randomly 
and indicate that the edge trajectory has entered the (phase-space) vicinity of an exact 
periodic solution, at which Tmin vanishes. The starting point at t ~ 150 which yielded the 
lowest rmin was used as an initial state for the Newton-Krylov algorithm with Az = 27r/a 
and A6' = 0. This readily converged to the asymmetric TW, shifted azimuthally by a half 
turn, with Tmin ~ 10~^°. 

Shooting in the opposite direction, then refining properly the angle 0, leads to another 
'edge' trajectory H2. The corresponding recurrence function rminif) is plotted in Figure 
[TOl The features of Tmin are reminiscent of those of the trajectory Hi. A dip (labelled N2) 
appears near t k, 190 d/U , and the corresponding state has been used as a starting point for 
the Newton-Krylov algorithm (again we look explicitly for a TW solution). The algorithm 
converged to another TW solution, that we call 26_1.5. Close examination revealed that 
it is equivalent of the TW 26_1.25 mentioned in Kerswell & Tutty^^ and Duguet et. al.— , 
but with an axial wavenumber of a = 1.5 instead of a = 1.25. Here a = 0.75 = 1.5/2, 
so that exactly two wavelengths of the TW fill the whole domain. This TW solution 
was already mentioned in Duguet et. al.— because it is an attractor for edge trajectories 
when Re = 2400, a = 1.25, and when the flow is constrained to a 2-fold rotationally 
symmetric subspace. This attracting property is lost when no rotational symmetry is 
forced, which explains why the TW here first attracts and then repels the edge trajectory H2. 

Given the difficulty of getting the Newton-Krylov algorithm to converge in high dimen- 
sions unless the starting guess is sufficiently close to a solution, it is reasonable to assume 
that the above trajectories enter the neighbourhood of the TWs before ultimately relami- 
narising. Since the unstable manifold of these TWs in the laminar-turbulent boundary is 
of such small dimension, it is tempting to speculate that these visits are actually indicative 

21 



0.001 




0.01 



FIG. 10: Recurrence function rmin vs. time for the two 'edge' trajectories Hi (red solid) and 
H2 (blue dotted), starting from the perturbed RPO and constrained to the laminar-turbulent 
boundary. The dip labelled Nl (resp.N2) near t ~ 150 (resp. t ^ 190) indicates an approach 
towards a TW state. 



of heteroclinic connections linking the RPO to the two TWs. Establishing this formally is 
a significant undertaking not pursued further here but the idea that trajectories link differ- 
ent coherent structures to produce a saddle structure on the laminar-turbulent boundary is 
consistent with other recent observations.— 

The trajectories Hi and H2 are projected on a two-dimensional (J, D) diagram in Figure 
[TTl I and D are respectively the energy input and the total viscous dissipation in the 
computational domain, defined by / := ^ pu ■ ndS, D := Re~^ J |V x u\'^d^x. These two 
quantities are related by the energy equation 



dE 
It 



D. 



(27) 



where E is the total kinetic energy E := ^ J \u\'^d'^x . Hence (relative) equilibria are all 
located on the diagonal I = D. I and D are normalised in Figure [TT] by /, which is the value 
of / for the laminar flow (hence located here at (/, -D) = (1,1)). Both trajectories Hi and H2 
escape from the RPO solution along its unstable manifold, in a direction locally tangent to 
the laminar-turbulent boundary. They reach much larger values of D and J, before turning 
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FIG. 11: Normalised (/ — D) projection of the 'edge' trajectories ii\ (red) and i/2 (blue) starting 
from the perturbed RPO, Re = 1783.6. Also shown are the two TW states (blue cross and dot) 
identified using the Newton-Krylov algorithm : TWA stands for the Asymmetric TW, TW2b for 
the TW solution 26_1.5. The location of the RPO is undistinguishable from that of the Asymmetric 
TW on this plot. 

back towards the laminar state, each of them passing by the vicinity of one of the TWs. 
The time of closest approach in the (J, D) plane does not coincide with the minima of r^m 
(e.g. N2 on H2 is not the closest point of approach to 26_1.5 in the {I,D) plane). This 
difference, of course, highlights the ongoing dilemma of how to select the 'right' norm to 
measure distances in the infinite-dimensional phase space associated to the Navier-Stokes 
equations. 



23 



V. CONCLUSION 

We have shown that it is possible to numerically capture relative periodic orbits (RPOs) 
in pipe flow when the associated dynamical system size exceeds 10^ degrees of freedom. 
We have applied successfully our method to find a RPO branch in the vicinity of the Hopf 
bifurcation occurring along the asymmetric Travelling Wave branch for a = 0.75. This 
exact coherent structure consists of a streaky pattern which is modulated over one period 
of roughly 43 d/U, while travelling approximately 60 pipe diameters downwards in the axial 
direction. 

Further evidence is presented here that the laminar-turbulent boundary is structured 
around of network of exact solutions linked to each other by heteroclinic and homoclinic 
connect ions^'i^i^Si'^. For L ^ 5d {a = 0.625) and Re = 2875, asymmetric TW solutions are 
recurrently approached by laminar-turbulent trajectories;— here for L ~ 4.19(i (a = 0.75) 
and Re = 1783.6, the same conclusion seems to hold. However, now this schematic view of 
edge trajectories recurrently visiting neighbourhoods of TWs needs to be extended to take 
into account relative periodic orbits. All these solutions and their stable manifolds make 
up the laminar-turbulent boundary. 

The main problem in accurately identifying RPOs is the exponential divergence of a 
nearby orbit from the RPO. In practice, the instability of the RPO together with its typically 
long period lead to a large value for its leading characteristic multiplier. This causes a 
significant loss of numerical accuracy when searching for neighbouring orbits which almost 
close on themselves. This effect was found to worsen away from the bifurcation point as the 
RPO becomes more unstable so that it was only possible to trace the RPO branch accurately 
in a relatively small neighbourhood of the bifurcation point. As a result, its amplitude 
remains small and the RPO dynamics is close to that of the underlying TW solution. The 
original goal of tracing the RPO down (in Re) to a likely saddle node bifurcation point 
and then exploring the 'other' branch have, unfortunately, been out of reach. Whether 
this 'other' branch reconnects to another TW branch or becomes a more nonlinear unstable 
object possibly embedded in the flow's turbulent dynamics remains a challenge for the future. 
There are ways to improve the numerical algorithm such as a multiple-shooting technique 
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where the orbit is subdivided into manageable pieces which now clearly need to be explored. 
Physically, the instability of the RPO means that phase-space trajectories are unlikely 
to spend a long time in its vicinity, despite the likelihood that the RPO is linked via 
heteroclinic connections to other more attracting solutions. In particular, the probability is 
low that such a solution could be identified experimentally, despite recent progress made in 
fiow visualisation or indeed numerically as part of a direct numerical simulation. The same 
conclusion may well apply to all RPOs of the system. 

At the start of this study, the stability of various lower TW branches was considered in 
order to identify how frequently Hopf bifurcations of TW branches arise. For a = 0.75, 
the asymmetric TW branch undergoes only one Hopf bifurcation below Re = 5000 which is 
examined in this paper. The mirror-symmetric branch displays another one near Re = 3487, 
its period near the bifurcating point is T = 68, which makes it even more unstable and harder 
to continue than the previous one. For the same parameters, we have checked that another 
known branch of TW (referred to as 26_1.25^'^) does not display any Hopf bifurcations at 
all below Re = 3400. Despite the fact that many more TWs might exist, it appears that 
Hopf bifurcations of lower-branch TW solutions are not frequent, and that short-period 
RPOs associated with them are not generic objects. From this point of view, the fact that 
mainly TW solutions have been identified during transition does not mean that RPOs do 
not participate in the transitional dynamics, but rather that their contribution is likely to be 
more infrequent and fieeting. However, it might not be the case for other shear fiows such as 
plane Couette fiow where (non-relative) periodic orbits also bifurcate from a lower-branch 
steady state.— "^"^ This does not exclude either the existence of RPO solutions entirely 
disconnected from the TW solutions, which would appear via saddle-node bifurcations. 
Such a (periodic) solution has been computed in plane Couette fiow, and its dynamics is 
embedded in the turbulent dynamics.— i^ The numerical method described in Section |TT] is 
an adequate tool to seek them providing a good enough starting point is identified for the 
Newton-Krylov algorithm. 
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